/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
    \file mesh.cpp
    \brief This file contains the Mesh class

   The Mesh class contains data and methods for the rappresentation
   and the use of a two dimensional SEM mash over a rectangular domain
*/

# include "mesh.hpp"

namespace spectral {

Mesh::Mesh(int nx,
           int ny,
           int nnodes,
           int nmax,
           Point<double> const & A,
           Point<double> const & B,
           double legDelta,
           double legEps,
           int order) : _nx(nx),
                        _ny(ny),
                        _nnodes(nnodes),
                        _A(A),
                        _B(B),
                        _legDelta(legDelta),
                        _legEps(legEps),
                        _integr(nnodes,legDelta,legEps,nmax,order) {
    assert(_nx>0);
    assert(_ny>0);
    assert(_legDelta>0);
    assert(_legEps>0);

    _vertexX=new double[_nx+1];
    _vertexY=new double[_ny+1];

    _vertexX[0]=_A.get(0);
    _vertexX[_nx]=_B.get(0);
    _vertexY[0]=_A.get(1);
    _vertexY[_ny]=_B.get(1);

    double stepX=(_vertexX[_nx]-_vertexX[0])/double(_nx);
    double stepY=(_vertexY[_ny]-_vertexY[0])/double(_ny);

    assert(stepX>0);
    assert(stepY>0);

    for(int i=1;i<_nx;++i) {
        _vertexX[i]=_vertexX[i-1]+stepX;
    }
    for(int i=1;i<_ny;++i) {
        _vertexY[i]=_vertexY[i-1]+stepY;
    }
}

Mesh::Mesh(Mesh const & Me): _nx(Me._nx),
                             _ny(Me._ny),
                             _nnodes(Me._nnodes),
                             _A(Me._A),
                             _B(Me._B),
                             _legDelta(Me._legDelta),
                             _legEps(Me._legEps),
                             _integr(Me._integr) {
    _vertexX=new double[_nx+1];
    _vertexY=new double[_ny+1];

    for(int i=0;i<_nx+1;++i) {
        _vertexX[i]=Me._vertexX[i];
    }
    for(int i=0;i<_ny+1;++i) {
        _vertexY[i]=Me._vertexY[i];
    }
}

Mesh::Mesh(Mesh const & Me,
           Point<int> const & A,
           Point<int> const & B): _nnodes(_nnodes),
                                  _legDelta(_legDelta),
                                  _legEps(_legEps),
                                  _integr(Me._integr) {
    int Ax=A.get(0);
    int Ay=A.get(1);
    int Bx=B.get(0);
    int By=B.get(1);

    assert(Ax < Bx && Ay < By);
    assert(Ax >= 0 && Bx <= Me._nx);
    assert(Ay >= 0 && By <= Me._ny);

    _nx=Bx-Ax;
    _ny=By-Ay;
    _vertexX=new double[_nx+1];
    _vertexY=new double[_ny+1];
    for(int i=0;i<_nx+1;++i) {
        _vertexX[i]=Me._vertexX[i+Ax];
    }
    for(int i=0;i<_ny+1;++i) {
        _vertexY[i]=Me._vertexY[i+Ay];
    }
}

Mesh::Mesh(Mesh const & Me,
           Point<double> const & A,
           Point<double> const & B): _nnodes(_nnodes),
                                     _legDelta(_legDelta),
                                     _legEps(_legEps),
                                     _integr(Me._integr) {
    Point<int> AI(gridPosition(A));
    Point<int> BI(gridPosition(B));
    Mesh(Me,AI,BI);
}

Mesh::~Mesh() {
    delete[] _vertexX;
    delete[] _vertexY;
}

double Mesh::integrate(std::valarray<double> const & v,
                       int elemnum) const {
    Point<double> C;
    Point<double> D;
    elementNumber(elemnum,C,D);

    double * val;
    std::vector<int> nodes;
    std::vector<int>::const_iterator it;

    getGlobalNum(elemnum,nodes);
    val = new double[_nnodes*_nnodes];

    int k(0);
    for(it=nodes.begin();it<nodes.end();++it) {
        val[k++]=v[*it];
    }

    double out=_integr.integrate(val,C,D);
    delete[] val;
    return out;
}

double Mesh::numNorm(std::valarray<double> const & v) const {
    double sum(0.0);
    std::valarray<double> tmp;

    tmp.resize(v.size());
    tmp=v*v;

    for(int i=0;i<_nx*_ny;++i) {
        sum+=integrate(tmp,i);
    }

    return std::sqrt(sum);
}

int Mesh::getnx() const{
    return _nx;
}

int Mesh::getny() const {
    return _ny;
}

int Mesh::getnelem() const {
    return _nx*_ny;
}

int Mesh::getnnodes( void ) const {
    return _nnodes;
}

int Mesh::getProblemDim() const {
    return (_nx*_nnodes-(_nx-1))*(_ny*_nnodes-(_ny-1));
}

int Mesh::getnbc() const {
    return 2*(_nx*_nnodes-(_nx-1))+2*(_ny*_nnodes-(_ny-1));
}

int Mesh::getnNodesPerEl() const {
    return _nnodes*_nnodes;
}

int Mesh::getnNodesPerBord(int border) const {
    switch(border) {
        case 0:
            return (_nx*_nnodes-(_nx-1));
        case 1:
            return (_ny*_nnodes-(_ny-1));
        case 2:
            return (_nx*_nnodes-(_nx-1));
        case 3:
            return (_ny*_nnodes-(_ny-1));
        default:
            exit(1);
    }
}

void Mesh::getElement(int elemnum,
                      Point<double> & BL,
                      Point<double> & TR) const {

    int nxcell=elemnum%_nx;
    int nycell=elemnum/_nx;

    assert( nxcell>=0 && nxcell<_nx );
    assert( nycell>=0 && nycell<_ny );

    BL[0]=_vertexX[nxcell];    BL[1]=_vertexY[nycell];
    TR[0]=_vertexX[nxcell+1];  TR[1]=_vertexY[nycell+1];
}

int Mesh::getGlobalNum(int elemnum,
                       int idxNode ) const {
    int nxcell=elemnum%_nx;
    int nycell=elemnum/_nx;
    int idx = idxNode%_nnodes;
    int idy = idxNode/_nnodes;

    // Nodes index at the beginning of the element
    int indx=_nnodes*nxcell-nxcell;
    int indy=_nnodes*nycell-nycell;

    int pdx=indx+idx;
    int pdy=indy+idy;

    int nodes_per_row=_nnodes*_nx-_nx+1;

    return pdx+pdy*nodes_per_row;
}

int Mesh::getGlobalNum(int elemnum,
                       int idx,
                       int idy) const {
    int nxcell=elemnum%_nx;
    int nycell=elemnum/_nx;

    // Nodes index at the beginning of the element
    int indx=_nnodes*nxcell-nxcell;
    int indy=_nnodes*nycell-nycell;

    int pdx=indx+idx;
    int pdy=indy+idy;

    int nodes_per_row=_nnodes*_nx-_nx+1;

    return pdx+pdy*nodes_per_row;
}

void Mesh::getGlobalNum(int elemnum,
                        std::vector<int> & nodes) const {
    assert( elemnum>=0 && elemnum<_nx*_ny );

    int numnodes=_nnodes*_nnodes;
    nodes.clear();
    nodes.reserve(numnodes);

    int nxcell=elemnum%_nx;
    int nycell=elemnum/_nx;

    // Nodes index at the beginning of the element
    int idx=_nnodes*nxcell-nxcell;
    int idy=_nnodes*nycell-nycell;

    // Nodes index at the end of the element+1
    int fdx=_nnodes*(nxcell+1)-nxcell;

    // Nodes per row
    int nodes_per_row=_nnodes*_nx-_nx+1;

    // Position
    int pdx(idx);
    int pdy(idy);

    for(int i=0;i<numnodes;++i) {
        nodes.push_back(pdx+pdy*nodes_per_row);
        pdx=( ( ++pdx<fdx || ++pdy<0 ) ? (pdx) : (idx) );
    }
}

int Mesh::getNodeNumber(int idx, int idy ) const {
    return idx+idy*getnNodesPerBord(0);
}

void Mesh::getBordElements(int border,
                           std::vector<int> & elements) const {
    assert( border>=0 && border<4 );

    elements.clear();

    switch(border) {
        case 0:
            elements.reserve(_nx);
            for(int i=0;i<_nx;++i) {
                elements.push_back(i);
            }
            break;
        case 1:
            elements.reserve(_ny);
            for(int i=0;i<_ny;++i) {
                elements.push_back(_nx*(i+1)-1);
            }
            break;
        case 2:
            elements.reserve(_nx);
            for(int i=0;i<_nx;++i) {
                elements.push_back(_nx*(_ny-1)+i);
            }
            break;
        case 3:
            elements.reserve(_ny);
            for(int i=0;i<_ny;++i) {
                elements.push_back(_nx*i);
            }
            break;
        default:
            exit(1);
    }
}

void Mesh::getBordNodes(int border,
                        std::vector<int> & nodes) const {

    assert( border>=0 && border<4 );

    nodes.clear();
    nodes.reserve(_nnodes);

    switch(border) {
        case 0:
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(i);
            }
            break;
        case 1:
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back((i+1)*_nnodes-1);
            }
            break;
        case 2:
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(_nnodes*(_nnodes-1)+i);
            }
            break;
        case 3:
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(i*_nnodes);
            }
            break;
        default:
            exit(1);
    }
}

void Mesh::getDomainBordNodes(int border,
			      std::vector<int> & nodes) const {
    nodes.clear();
    Point<int> P1(getVertexPosition(0, 0));
    Point<int> P2(getVertexPosition(_nx, _ny));
    
    switch(border) {
        case 0:
	    nodes.reserve(P2[0]-P1[0]+1);
	    for (int i(P1[0]); i<=P2[0]; ++i) {
	        nodes.push_back(getNodeNumber(i, P1[1]));
	    }
	    break;
	
	case 1:
	    nodes.reserve(P2[1]-P1[1]+1);
	    for (int i(P1[1]); i<=P2[1]; ++i) {
	        nodes.push_back(getNodeNumber(P2[0], i));
	    }
	    break;

	case 2:
	    nodes.reserve(P2[0]-P1[0]+1);
	    for (int i(P2[0]); i>=P1[0]; --i) {
	        nodes.push_back(getNodeNumber(i, P2[1]));
	    }
	    break;

	case 3:
	    nodes.reserve(P2[1]-P1[1]+1);
	    for (int i(P2[1]); i>=P1[1]; --i) {
	        nodes.push_back(getNodeNumber(P1[0], i));
	    }
	    break;

	default:
	    exit(1);
    }
}

void Mesh::whichBord(int node,
                     bool isonborder[4]) const {
    Point<int> tmp(gridPosition(node));
    for(int i=0;i<4;++i) {
        isonborder[i]=false;
    }
    if(tmp[1]==0) {
        isonborder[0]=true;
    }
    if(tmp[0]==_nx*(_nnodes-1)) {
        isonborder[1]=true;
    }
    if(tmp[1]==_ny*(_nnodes-1)) {
        isonborder[2]=true;
    }
    if(tmp[0]==0) {
        isonborder[3]=true;
    }
}

bool Mesh::isPhysicalVertex(int node) const {
    // maximum node number (not the number of nodes!)
    int nnx=_nx*(_nnodes-1);
    int nny=_ny*(_nnodes-1);
    Point<int> tmp(gridPosition(node));
    return ( tmp[0]==0   && tmp[1]==0   ||
             tmp[0]==nnx && tmp[1]==0   ||
             tmp[0]==nnx && tmp[1]==nny ||
             tmp[0]==0   && tmp[1]==nny );
}

void Mesh::getPhysicalVertex(int node,int & bord1,int & bord2) const {
    assert(isPhysicalVertex(node));
    bool isonborder[4];
    whichBord(node,isonborder);
    bord1=0;
    bord2=0;
    while(!isonborder[bord1]) {
        ++bord1;
    }
    bord2=bord1+1;
    while(!isonborder[bord2]) {
        ++bord2;
    }
    assert(bord1>=0 && bord1<4);
    assert(bord2>=0 && bord2<4);
}

int Mesh::getBordNum(int border,int i) const {
    switch(border) {
        case 0:
            return i;
        case 1:
            return i/(_nx*_nnodes-(_nx-1));
        case 2:
            return (_nx*_nnodes-(_nx-1))-1-i%(_nx*_nnodes-(_nx-1));
        case 3:
            return (_ny*_nnodes-(_ny-1))-1-i/(_nx*_nnodes-(_nx-1));
        default:
            exit(1);
    }
}

void Mesh::getGlobBordNodes(int border,
                            int elemnum,
                            std::vector<int> & nodes) const {
    assert( border>=0 && border<4 );

    nodes.clear();
    nodes.reserve(_nnodes);

    // Element position
    int nxcell=elemnum%_nx;
    int nycell=elemnum/_nx;

    // Nodes index at the beginning of the element
    int idx=_nnodes*nxcell-nxcell;
    int idy=_nnodes*nycell-nycell;

    // Nodes index at the end of the element
    int fdx=_nnodes*(nxcell+1)-nxcell-1;
    int fdy=_nnodes*(nycell+1)-nycell-1;

    // Nodes per row
    int nodes_per_row=_nnodes*_nx-_nx+1;

    int index;

    switch(border) {
        case 0:
            index=idx;
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(index+idy*nodes_per_row);
                ++index;
            }
            break;
        case 1:
            index=idy;
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(fdx+index*nodes_per_row);
                ++index;
            }
            break;
        case 2:
            index=idx;
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(index+fdy*nodes_per_row);
                ++index;
            }
            break;
        case 3:
            index=idy;
            for(int i=0;i<_nnodes;++i) {
                nodes.push_back(idx+index*nodes_per_row);
                ++index;
            }
            break;
        default:
            exit(1);
    }
}

Point<double> Mesh::getNodePosition(int node) const {
    Point<int> tmp(gridPosition(node));
    return getNodePosition(tmp[0],tmp[1]);
}

Point<double> Mesh::getNodePosition(int pdx,int pdy) const {
    // Element number (we give the last node to the next element
    // except for the last element of a row/column)
    Point<double> tmp=getRefNodePosition(pdx,pdy);

    int edx=std::min(pdx/(_nnodes-1),_nx-1);
    int edy=std::min(pdy/(_nnodes-1),_ny-1);

    Point<double> OUT;
    Point<double> BL;
    Point<double> TR;

    BL[0]=_vertexX[edx];    BL[1]=_vertexY[edy];
    TR[0]=_vertexX[edx+1];  TR[1]=_vertexY[edy+1];

    OUT[0]=0.5*(TR[0]-BL[0])*(tmp[0]+1.0)+BL[0];
    OUT[1]=0.5*(TR[1]-BL[1])*(tmp[1]+1.0)+BL[1];

    return OUT;
}

Point<double> Mesh::getRefNodePosition(int node) const {
    Point<int> tmp(gridPosition(node));
    return getRefNodePosition(tmp[0],tmp[1]);
}

Point<double> Mesh::getRefNodePosition(int pdx,int pdy) const {
    int edx=std::min(pdx/(_nnodes-1),_nx-1);
    int edy=std::min(pdy/(_nnodes-1),_ny-1);

    // Position in the local coordinates
    int idx=pdx-edx*_nnodes+edx;
    int idy=pdy-edy*_nnodes+edy;

    Point<double> OUT;
    OUT[0]=_integr.getnode(idx);
    OUT[1]=_integr.getnode(idy);
    return OUT;
}

Point<int> Mesh::getVertexPosition(int edx,int edy) const {
    Point<int> P;
    P[0]=edx*(_nnodes-1);
    P[1]=edy*(_nnodes-1);
    return P;
}

Point<int> Mesh::gridPosition(int node) const {
    // Number of nodes along the x-axis
    int nnx=_nx*_nnodes-_nx+1;

    // Position in the nodes grid
    int pdx=node%nnx;
    int pdy=node/nnx;

    return Point<int>(pdx,pdy);
}

// naive implementation (alpha version) (this is almost useless)
Point<int> Mesh::gridPosition(Point<double> const & P) const {
    Point<int> I;

    double Px=P.get(0);
    double Py=P.get(1);

    double deltaold=1e16;
    double delta=std::abs(Px-_vertexX[0]);

    int i=1;

    while(delta<deltaold && i<_nx+1) {
        deltaold=delta;
        delta=std::abs(Px-_vertexX[i]);
        ++i;
    }

    I[0]=i-1;

    i=1;
    deltaold=1e16;
    delta=std::abs(Py-_vertexY[0]);

    while(delta<deltaold && i<_ny+1) {
        deltaold=delta;
        delta=std::abs(Py-_vertexY[i]);
        ++i;
    }

    I[1]=i-1;

    return I;
}


int Mesh::whereami(Point<double> const & P) const {
    Point<int> I;

    double Px=P.get(0);
    double Py=P.get(1);

    assert(Px>=_A.get(0) && Px<=_B.get(0));
    assert(Py>=_A.get(1) && Py<=_B.get(1));

    double stepX=(_B.get(0)-_A.get(0))/double(_nx);
    double stepY=(_B.get(1)-_A.get(1))/double(_ny);

    int idx=(int)(Px/stepX);
    int idy=(int)(Py/stepY);

    return idx+idy*_nx;
}

void Mesh::elementNumber(int n,
                         Point<double> & P1,
                         Point<double> & P2) const {

    int nxcell=n%_nx;
    int nycell=n/_nx;

    assert(n >= 0 && n < _nx*_ny);
    assert(nxcell < _nx);
    assert(nycell < _ny);

    P1[0]=_vertexX[nxcell];
    P1[1]=_vertexY[nycell];
    P2[0]=_vertexX[nxcell+1];
    P2[1]=_vertexY[nycell+1];
}

void Mesh::save_to_file(std::valarray<double> const & vec,
                        char const * filename) const {
    save_to_file(&vec[0],filename);
}

void Mesh::save_to_file(double const * vec,
                        char const * filename) const {
    std::ofstream outfile;
    outfile.open(filename,std::ios_base::out);

    Point<double> P;
    int _vertex;

    // Number of nodes along the x-axis
    int nnx=_nx*_nnodes-_nx+1;
    // Number of nodes along the y-axis
    int nny=_ny*_nnodes-_ny+1;

    for(int idx=0;idx<nnx;++idx) {
        for(int idy=0;idy<nny;++idy) {
             _vertex=idx+idy*nnx;
             P=getNodePosition(idx,idy);
             outfile << P[0] << " " << P[1] << " " << vec[_vertex] << std::endl;
        }
        outfile << std::endl;
    }
    outfile.close();
}

} // namespace
